module riemann_mod

  use spark_cubed_sphere
  use flogger
  use params_mod
  use physics_const_mod

  implicit none

  private

  public riemann_init
  public riemann_solver_adv
  public riemann_solver_sw
  public riemann_solver_hor
  public riemann_solver_ver
  public calc_pressure

  interface
    function riemann_solver_adv_interface(u, qL, qR, ghs)
      import r_kind
      real(r_kind), intent(in) :: u
      real(r_kind), intent(in) :: qL(:)
      real(r_kind), intent(in) :: qR(:)
      real(r_kind), intent(in) :: ghs
      real(r_kind) riemann_solver_adv_interface(size(qL))
    end function riemann_solver_adv_interface
    
    function riemann_solver_sw_interface(axis, iG, sqrtG, qL, qR, ghsL, ghsR)
      import r_kind
      integer , intent(in) :: axis
      real(r_kind), intent(in) :: iG(2,2)
      real(r_kind), intent(in) :: sqrtG
      real(r_kind), intent(in) :: qL(:)
      real(r_kind), intent(in) :: qR(:)
      real(r_kind), intent(in) :: ghsL
      real(r_kind), intent(in) :: ghsR
      real(r_kind) riemann_solver_sw_interface(size(qL))
    end function riemann_solver_sw_interface
    
    function riemann_solver_hor_interface(axis, sqrtGhL, sqrtGhR, iGhL, iGhR, qL, qR)
      import r_kind
      integer , intent(in) :: axis
      real(r_kind), intent(in) :: sqrtGhL
      real(r_kind), intent(in) :: sqrtGhR
      real(r_kind), intent(in) :: iGhL(2,2)
      real(r_kind), intent(in) :: iGhR(2,2)
      real(r_kind), intent(in) :: qL(:)
      real(r_kind), intent(in) :: qR(:)
      real(r_kind) riemann_solver_hor_interface(size(qL))
    end function riemann_solver_hor_interface
    
    function riemann_solver_ver_interface(sqrtGL, sqrtGR, sqrtGhL, sqrtGhR, sqrtGvL, sqrtGvR, iGL, iGR, iGvL, iGvR, qL, qR)
      import r_kind
      real(r_kind), intent(in) :: sqrtGL
      real(r_kind), intent(in) :: sqrtGR
      real(r_kind), intent(in) :: sqrtGhL
      real(r_kind), intent(in) :: sqrtGhR
      real(r_kind), intent(in) :: sqrtGvL
      real(r_kind), intent(in) :: sqrtGvR
      real(r_kind), intent(in) :: iGL(3,3)
      real(r_kind), intent(in) :: iGR(3,3)
      real(r_kind), intent(in) :: iGvL(2,3:3)
      real(r_kind), intent(in) :: iGvR(2,3:3)
      real(r_kind), intent(in) :: qL(:)
      real(r_kind), intent(in) :: qR(:)
      real(r_kind) riemann_solver_ver_interface(size(qL))
    end function riemann_solver_ver_interface
  end interface

  procedure(riemann_solver_adv_interface), pointer :: riemann_solver_adv=> null()
  procedure(riemann_solver_sw_interface ), pointer :: riemann_solver_sw => null()
  procedure(riemann_solver_hor_interface), pointer :: riemann_solver_hor=> null()
  procedure(riemann_solver_ver_interface), pointer :: riemann_solver_ver=> null()

contains

  subroutine riemann_init()

    select case(model_type)
    case(0)
      riemann_solver_adv => riemann_solver_adv_upwind
    case(1)
      select case (riemann_scheme_hor)
        case ('AUSM_up')
          riemann_solver_sw => riemann_solver_sw_ausm_up
        case ('LLF')
          riemann_solver_sw => riemann_solver_sw_LLF
        case default
          call log_error('Unknown riemann_scheme_hor "' // trim(riemann_scheme_hor) // '"!', pid=proc%id)
      end select
    case(2)
      select case (riemann_scheme_hor)
        case ('AUSM_up')
          riemann_solver_hor => riemann_solver_hor_ausm_up
        !case ('LLF')
        !  riemann_solver_hor => riemann_solver_hor_LLF
        case default
          call log_error('Unknown riemann_scheme_hor "' // trim(riemann_scheme_hor) // '"!', pid=proc%id)
        end select
        
      select case (riemann_scheme_ver)
        case ('AUSM_up')
          riemann_solver_ver => riemann_solver_ver_ausm_up
        !case ('LLF')
        !  riemann_solver_ver => riemann_solver_ver_LLF
        case default
          call log_error('Unknown riemann_scheme_ver "' // trim(riemann_scheme_ver) // '"!', pid=proc%id)
      end select
    end select
    
  end subroutine riemann_init

  function riemann_solver_adv_upwind(u, qL, qR, ghs) result(res)
    real(r_kind), intent(in) :: u
    real(r_kind), intent(in) :: qL(:)
    real(r_kind), intent(in) :: qR(:)
    real(r_kind), intent(in) :: ghs
    real(r_kind) res(size(qL))
    
    res = u * 0.5_r_kind * (qL + qR - sign(1.0_r_kind, u) * (qR - qL))
    
  end function riemann_solver_adv_upwind

  function riemann_solver_sw_ausm_up(axis, iG, sqrtG, qL, qR, ghsL, ghsR) result(res)

    integer , intent(in) :: axis
    real(r_kind), intent(in) :: iG(2,2)
    real(r_kind), intent(in) :: sqrtG
    real(r_kind), intent(in) :: qL(:)
    real(r_kind), intent(in) :: qR(:)
    real(r_kind), intent(in) :: ghsL
    real(r_kind), intent(in) :: ghsR
    real(r_kind) res(size(qL))

    real(r_kind) ghL, ghR, gh, sqrtIG
    real(r_kind) pL, pR, p         ! Pressure
    real(r_kind) uL, uR            ! Flow speed
    real(r_kind) cL, cR, c         ! Wave speed
    real(r_kind) ML, MR, Mb2, Ma   ! Mach number
    real(r_kind) m, p5p, p5n

    real(r_kind), parameter :: Ku  = 0.75_r_kind
    real(r_kind), parameter :: Kp  = 0.25_r_kind
    real(r_kind), parameter :: sigma =  1.0_r_kind

    sqrtIG = sqrt(iG(axis,axis))

    ghL = qL(1) / sqrtG + ghsL
    ghR = qR(1) / sqrtG + ghsR
    gh  = 0.5_r_kind * (ghL + ghR)

    pL  = 0.5_r_kind * ghL**2
    pR  = 0.5_r_kind * ghR**2

    uL  = qL(1+axis) / qL(1) / sqrtIG
    uR  = qR(1+axis) / qR(1) / sqrtIG

    cL  = sqrt(ghL)
    cR  = sqrt(ghR)
    c   = 0.5_r_kind * (cL + cR)

    ML  = uL / c
    MR  = uR / c

    Mb2 = (uL**2 + uR**2) / (2 * c**2)

    Ma = M4(ML, 1) + M4(MR, -1) - Kp * max(1.0_r_kind - sigma * Mb2, 0.0_r_kind) * (ghR**2 - ghL**2) / (2 * gh * c**2)
    m = c * Ma * sqrtIG

    p5p = P5(ML,  1)
    p5n = P5(MR, -1)
    p = p5p * pL + p5n * pR - Ku * p5p * p5n * (ghL + ghR) * c * (uR - uL)

    res = m * 0.5_r_kind * (qL + qR - sign(1.0_r_kind, m) * (qR - qL))
    res(2) = res(2) + sqrtG * iG(1,axis) * p
    res(3) = res(3) + sqrtG * iG(2,axis) * p
    
  end function riemann_solver_sw_ausm_up
  
  function riemann_solver_nh(rhoL, rhoR, pL, pR, upL, upR, u1L, u1R, u2L, u2R, ptL, ptR, mrL, mrR) result(res)
    real(r_kind), intent(in) :: rhoL
    real(r_kind), intent(in) :: rhoR
    real(r_kind), intent(in) :: pL
    real(r_kind), intent(in) :: pR
    real(r_kind), intent(in) :: upL
    real(r_kind), intent(in) :: upR
    real(r_kind), intent(in) :: u1L
    real(r_kind), intent(in) :: u1R
    real(r_kind), intent(in) :: u2L
    real(r_kind), intent(in) :: u2R
    real(r_kind), intent(in) :: ptL
    real(r_kind), intent(in) :: ptR
    real(r_kind), intent(in) :: mrL
    real(r_kind), intent(in) :: mrR
    real(r_kind) res(nVar)

    real(r_kind) rho             ! Density
    real(r_kind) p(nVar)         ! Pressure
    real(r_kind) cL, cR, a       ! Wave speed
    real(r_kind) ML, MR, Mb2, Ma ! Mach number
    real(r_kind) m, p5p, p5n

    real(r_kind) qL(nVar)
    real(r_kind) qR(nVar)
    
    real(r_kind), parameter :: Ku    = 0.75_r_kind
    real(r_kind), parameter :: Kp    = 0.25_r_kind
    real(r_kind), parameter :: sigma =  1.0_r_kind
    
    integer iVar
    
    rho = 0.5 * ( rhoL + rhoR )
    
    cL  = calc_sound_speed_nh(rhoL, pL, mrL)
    cR  = calc_sound_speed_nh(rhoR, pR, mrR)
    a   = 0.5_r_kind * (cL + cR)

    ML  = upL / a
    MR  = upR / a

    Mb2 = ( upL**2 + upR**2 ) / ( 2 * a**2 )

    Ma = M4(ML, 1) + M4(MR, -1) - Kp * max( 1.0_r_kind - sigma * Mb2, 0.0_r_kind ) * ( pR - pL ) / ( rho * a**2 )
    m  = a * Ma
    
    p5p  = P5(ML,  1)
    p5n  = P5(MR, -1)
    p    = 0
    p(2) = p5p * pL + p5n * pR - Ku * p5p * p5n * 2 * rho * a * ( upR - upL )
    
    qL(1) = rhoL
    qL(2) = rhoL * upL
    qL(3) = rhoL * u1L
    qL(4) = rhoL * u2L
    qL(5) = rhoL * ptL
    qL(6) = rhoL * mrL
    
    qR(1) = rhoR
    qR(2) = rhoR * upR
    qR(3) = rhoR * u1R
    qR(4) = rhoR * u2R
    qR(5) = rhoR * ptR
    qR(6) = rhoR * mrR
    
    res = m * 0.5_r_kind * ( qL + qR - sign(1.0_r_kind, m) * ( qR - qL ) ) + p
    
  end function riemann_solver_nh

  function riemann_solver_hor_ausm_up(axis, sqrtGL, sqrtGR, iGhL, iGhR, qL, qR) result(res)

    integer , intent(in) :: axis
    real(r_kind), intent(in) :: sqrtGL
    real(r_kind), intent(in) :: sqrtGR
    real(r_kind), intent(in) :: iGhL(2,2)
    real(r_kind), intent(in) :: iGhR(2,2)
    real(r_kind), intent(in) :: qL(:)
    real(r_kind), intent(in) :: qR(:)
    real(r_kind) res(size(qL))

    real(r_kind) sqrtIGL
    real(r_kind) sqrtIGR
    real(r_kind) sqrtGrhoL
    real(r_kind) sqrtGrhoR
    real(r_kind) rhoL, rhoR, rho   ! Density
    real(r_kind) pL, pR, p(nVar)   ! Pressure
    real(r_kind) uL, uR            ! Flow speed
    real(r_kind) mrL, mrR          ! mixing ratio
    real(r_kind) cL, cR, a         ! Wave speed
    real(r_kind) ML, MR, Mb2, Ma   ! Mach number
    real(r_kind) m, p5p, p5n
    
    real(r_kind) :: coefL(2:3)
    real(r_kind) :: coefR(2:3)

    real(r_kind), parameter :: Ku  = 0.75_r_kind
    real(r_kind), parameter :: Kp  = 0.25_r_kind
    real(r_kind), parameter :: sigma =  1.0_r_kind
    
    integer iVar

    sqrtIGL = sqrt(iGhL(axis,axis))
    sqrtIGR = sqrt(iGhR(axis,axis))
    
    sqrtGrhoL = qL(1) + qL(6)
    sqrtGrhoR = qR(1) + qR(6)
    
    rhoL  = sqrtGrhoL / sqrtGL
    rhoR  = sqrtGrhoR / sqrtGR
    rho   = 0.5 * ( rhoL + rhoR )
    
    pL  = calc_pressure(sqrtGL,qL)
    pR  = calc_pressure(sqrtGR,qR)
    
    ! Convert wind from computational space to physical space by divide sqrtIG
    uL  = qL(1+axis) / sqrtGrhoL / sqrtIGL
    uR  = qR(1+axis) / sqrtGrhoR / sqrtIGR
    
    mrL = qL(6) / qL(1)
    mrR = qR(6) / qR(1)
    
    cL  = calc_sound_speed_nh(rhoL, pL, mrL)
    cR  = calc_sound_speed_nh(rhoR, pR, mrR)
    a   = 0.5_r_kind * (cL + cR)

    ML  = uL / a
    MR  = uR / a

    Mb2 = ( uL**2 + uR**2 ) / ( 2 * a**2 )

    Ma = M4(ML, 1) + M4(MR, -1) - Kp * max(1.0_r_kind - sigma * Mb2, 0.0_r_kind) * (pR - pL) / ( rho * a**2)
    m = 0.5 * ( cL * sqrtIGL + cR * sqrtIGR ) * Ma

    coefL(2) = sqrtGL * iGhL(1,axis)
    coefR(2) = sqrtGR * iGhR(1,axis)
    coefL(3) = sqrtGL * iGhL(2,axis)
    coefR(3) = sqrtGR * iGhR(2,axis)
    
    p5p = P5(ML,  1)
    p5n = P5(MR, -1)
    p   = 0
    do iVar = 2, 3
      p(iVar) = p5p * coefL(iVar) * pL + p5n * coefR(iVar) * pR - Ku * p5p * p5n * ( coefL(iVar) * rhoL + coefR(iVar) * rhoR ) * a * (uR - uL)
    enddo

    res = m * 0.5_r_kind * (qL + qR - sign(1.0_r_kind, m) * (qR - qL)) + p
    
  end function riemann_solver_hor_ausm_up

  function riemann_solver_ver_ausm_up(sqrtGL, sqrtGR, sqrtGhL, sqrtGhR, sqrtGvL, sqrtGvR, iGL, iGR, iGvL, iGvR, qL, qR) result(res)
    real(r_kind), intent(in) :: sqrtGL
    real(r_kind), intent(in) :: sqrtGR
    real(r_kind), intent(in) :: sqrtGhL
    real(r_kind), intent(in) :: sqrtGhR
    real(r_kind), intent(in) :: sqrtGvL
    real(r_kind), intent(in) :: sqrtGvR
    real(r_kind), intent(in) :: iGL(3,3)
    real(r_kind), intent(in) :: iGR(3,3)
    real(r_kind), intent(in) :: iGvL(2,3:3)
    real(r_kind), intent(in) :: iGvR(2,3:3)
    real(r_kind), intent(in) :: qL(:)
    real(r_kind), intent(in) :: qR(:)
    real(r_kind) res(size(qL))

    real(r_kind) sqrtGrhoL
    real(r_kind) sqrtGrhoR
    real(r_kind) rhoL, rhoR, rho   ! Density
    real(r_kind) pL, pR, p(nVar)   ! Pressure
    real(r_kind) uL, uR            ! Flow speed
    real(r_kind) mrL, mrR          ! mixing ratio
    real(r_kind) cL, cR, a         ! Wave speed
    real(r_kind) ML, MR, Mb2, Ma   ! Mach number
    real(r_kind) m, p5p, p5n
    
    real(r_kind) :: coefL(2:4)
    real(r_kind) :: coefR(2:4)

    real(r_kind), parameter :: Ku  = 0.75_r_kind
    real(r_kind), parameter :: Kp  = 0.25_r_kind
    real(r_kind), parameter :: sigma =  1.0_r_kind
    
    integer iVar

    sqrtGrhoL = qL(1) + qL(6)
    sqrtGrhoR = qR(1) + qR(6)
    
    rhoL = sqrtGrhoL / sqrtGL
    rhoR = sqrtGrhoR / sqrtGR
    rho  = 0.5 * ( rhoL + rhoR )
    
    pL  = calc_pressure(sqrtGL,qL)
    pR  = calc_pressure(sqrtGR,qR)
    
    uL  = ( qL(4) + sqrtGvL * ( iGvL(1,3) * qL(2) + iGvL(2,3) * qL(3) ) ) / sqrtGrhoL ! a sqrtGv has been multiplyed here
    uR  = ( qR(4) + sqrtGvR * ( iGvR(1,3) * qR(2) + iGvR(2,3) * qR(3) ) ) / sqrtGrhoR ! a sqrtGv has been multiplyed here

    mrL = qL(6) / qL(1)
    mrR = qR(6) / qR(1)
    
    cL  = calc_sound_speed_nh(rhoL, pL, mrL)
    cR  = calc_sound_speed_nh(rhoR, pR, mrR)
    a   = 0.5_r_kind * (cL + cR)
    
    ML  = uL / a
    MR  = uR / a

    Mb2 = ( uL**2 + uR**2 ) / ( 2 * a**2 )
    
    
    Ma = M4(ML, 1) + M4(MR, -1) - Kp * max(1.0_r_kind - sigma * Mb2, 0.0_r_kind) * (pR - pL) / ( rho * a**2 )
    m = 0.5 * ( cL / sqrtGvL + cR / sqrtGvR ) * Ma
    
    coefL(2) = sqrtGL * iGL(1,3)
    coefR(2) = sqrtGR * iGR(1,3)
    coefL(3) = sqrtGL * iGL(2,3)
    coefR(3) = sqrtGR * iGR(2,3)
    coefL(4) = sqrtGhL
    coefR(4) = sqrtGhR
    
    p5p = P5(ML,  1)
    p5n = P5(MR, -1)
    p   = 0
    do iVar = 2, 4
      p(iVar) = p5p * coefL(iVar) * pL + p5n * coefR(iVar) * pR - Ku * p5p * p5n * ( coefL(iVar) * rhoL + coefR(iVar) * rhoR ) * a * (uR - uL)
      !p(iVar) = coefL(iVar) * ( p5p * pL + p5n * pR - Ku * p5p * p5n * ( rhoL + rhoR ) * a * (uR - uL) )
    enddo
    
    res = m * 0.5_r_kind * (qL + qR - sign(1.0_r_kind, m) * (qR - qL)) + p
    
  end function riemann_solver_ver_ausm_up

  pure real(r_kind) function M2(M, s)
  
    real(r_kind), intent(in) :: M
    integer , intent(in) :: s
  
    M2 = s * 0.25_r_kind * (M + s)**2
  
  end function M2
  
  pure real(r_kind) function M4(M, signal)
  
    real(r_kind), intent(in) :: M
    integer     , intent(in) :: signal
    
    real(r_kind), parameter :: beta = 0.125_r_kind
    
    real(r_kind) :: M41, M42
    real(r_kind) :: s
    
    !if (abs(M) >= 1) then
    !  M4 = 0.5_r_kind * (M + signal * abs(M))
    !else
    !  M4 = M2(M, signal) * (1._r_kind - signal * 16 * beta * M2(M, -signal))
    !end if

    M41 = 0.5_r_kind * ( M + signal * abs(M) )
    M42 = M2( M, signal ) * ( 1._r_kind - signal * 16._r_kind * beta * M2( M, - signal ) )
    
    s = sig(abs(M),1._r_kind)
    
    M4 = s * M41 + ( 1._r_kind - s) * M42
  
  end function M4
  
  pure real(r_kind) function P5(M, signal)
  
    real(r_kind), intent(in) :: M
    integer , intent(in) :: signal
  
    real(r_kind), parameter :: alpha = 0.1875_r_kind
    real(r_kind) :: P51, P52
    real(r_kind) :: s, sm
  
    !if (abs(M) >= 1) then
    !  P5 = 0.5_r_kind * (1._r_kind + signal * sign(1.0_r_kind, M))
    !else
    !  P5 = M2(M, signal) * ((2._r_kind * signal - M) - signal * 16._r_kind * alpha * M * M2(M, -signal))
    !end if

    s = sig( M, 0._r_kind )
    sm = s - ( 1._r_kind - s )
    
    P51 = 0.5_r_kind * ( 1._r_kind + signal * sm )
    P52 = M2( M, signal ) * ( ( 2._r_kind * signal - M ) - signal * 16._r_kind * alpha * M * M2( M, -signal ) )
    
    s = sig( abs(M), 1._r_kind )
    
    P5 = s * P51 + ( 1._r_kind - s ) * P52
  
  end function P5
  
  pure real(r_kind) function sig(x,y)
    real(r_kind), intent(in) :: x, y
    real(r_kind), parameter :: K = tolerance
    
    sig = 0.5_r_kind * ( -( 1._r_kind + K ) * ( x - y ) / ( -abs( x - y ) - K ) + 1._r_kind );
    
  end function sig
  
  function riemann_solver_sw_LLF(axis, iG, sqrtG, qL, qR, ghsL, ghsR) result(res)
    integer , intent(in) :: axis
    real(r_kind), intent(in) :: iG(2,2)
    real(r_kind), intent(in) :: sqrtG
    real(r_kind), intent(in) :: qL(:)
    real(r_kind), intent(in) :: qR(:)
    real(r_kind), intent(in) :: ghsL
    real(r_kind), intent(in) :: ghsR
    real(r_kind) res(size(qL))
    
    real(r_kind) FL(size(qL))
    real(r_kind) FR(size(qR))
    real(r_kind) lambdaL, lambdaR
    real(r_kind) ghL, ghR
    real(r_kind) flowL, flowR
    
    FL = calc_F(axis,sqrtG,iG,qL,ghsL,ghL,flowL)
    FR = calc_F(axis,sqrtG,iG,qR,ghsR,ghR,flowR)
    lambdaL = calc_eignvalue(axis,iG,ghL,flowL)
    lambdaR = calc_eignvalue(axis,iG,ghR,flowR)
    
    res = 0.5 * ( FL + FR - max(lambdaL,lambdaR) * ( qR - qL ) )
    
  contains
    function calc_F(axis,sqrtG,iG,q,ghs,gh,flow) result(res)
      integer , intent(in) :: axis
      real(r_kind), intent(in ) :: sqrtG
      real(r_kind), intent(in ) :: iG(2,2)
      real(r_kind), intent(in ) :: q(:)
      real(r_kind), intent(in ) :: ghs
      real(r_kind), intent(out) :: gh
      real(r_kind), intent(out) :: flow
      real(r_kind) :: res(size(q))
      
      real(r_kind) :: u(2), p
      integer :: axis1
      
      axis1 = axis + 1
      gh   = q(1) / sqrtG + ghs
      u(1) = q(2) / q(1)
      u(2) = q(3) / q(1)
      p    = 0.5 * gh * gh
      flow = u(axis)
      
      res(1) = q(axis1)
      res(2) = q(axis1) * u(1) + sqrtG * iG(1,axis) * p
      res(3) = q(axis1) * u(2) + sqrtG * iG(2,axis) * p
      
    end function calc_F
    
    function calc_eignvalue(axis,iG,gh,flow) result(res)
      integer , intent(in) :: axis
      real(r_kind), intent(in) :: iG(2,2)
      real(r_kind), intent(in) :: gh
      real(r_kind), intent(in) :: flow
      real(r_kind) :: res
      
      res = abs( flow ) + sqrt( iG(axis,axis) * gh )
      
    end function calc_eignvalue
  end function riemann_solver_sw_LLF
  
  function calc_pressure(sqrtG,q)
    real(r_kind) calc_pressure
    real(r_kind) sqrtG
    real(r_kind) q(6)
    
    real(r_kind) w1
    real(r_kind) w5
    real(r_kind) w6
    
    real(r_kind) rho
    real(r_kind) R
    real(r_kind) gamma
    real(r_kind) theta
    real(r_kind) cp
    real(r_kind) cv
    real(r_kind) sh
    real(r_kind) kappa
    real(r_kind) sqrtGrho
    
    w1 = q(1)
    w5 = q(5)
    w6 = q(6)
    
    sqrtGrho = w1 + w6
    
    rho   = sqrtGrho / sqrtG
    gamma = w6 / w1
    sh    = gamma / ( 1. + gamma )
    R     = ( 1. + eqv * sh ) * Rd
    theta = w5 / w1
    cp    = cpd + ( cpv - cpd ) * sh
    cv    = cvd + ( cvv - cvd ) * sh
    kappa = cp / cv
    
    calc_pressure = p0 * ( rho * theta * R / p0 )**kappa
    
  end function calc_pressure
  
  function calc_sound_speed_nh(rho, p, mr) result(c)
    real(r_kind) :: c
    real(r_kind) :: rho
    real(r_kind) :: p
    real(r_kind) :: mr
    
    real(r_kind) :: kappa
    
    kappa = ( Cpd + mr * Cpv ) / ( Cvd + mr * Cvv )
    
    c = sqrt( kappa * p / rho )
    
  end function calc_sound_speed_nh
  
end module riemann_mod
